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Abstract 

In this work the issue of whether key energetic properties (nonlinear, exponential- 
type dissipation in the absence of forcing and long-term stability under conditions of 
time dependent loading) are automatically inherited by deforming reference descrip- 
tions is resolved. These properties are intrinsic to real flows and the conventional 
Navier-Stokes equations. A completely general reference description of an incom- 
pressible, Newtonian fluid, which reconciles the differences between opposing schools 
of thought in the literature is derived for the purposes of this investigation. 

The work subsequently focusses on establishing a class of time discretisations which 
inherit these self-same energetic properties, irrespective of the time increment em- 
ployed. The findings of this analysis have profound consequences for the use of 
certain classes of finite difference schemes in the context of deforming references. 
It is significant that many algorithms presently in use do not automatically inherit 
the fundamental qualitative features of the dynamics. An "updated" approach as a 
means of avoiding ever burgeoning deformation gradients and a still further simpli- 
fied implementation are further topics explored. 

Keywords: Energy Conservation; Incompressible, Newtonian Fluid; Completely General 
Reference Description; Arbitrary Lagrangian Eulerian; A.L.E.; Rigid Body in a Fluid; 
Free Surface; Finite Elements; New Poincare Inequality. 



1 Introduction 

Descriptions of fluid motion are conventionally based on the principles of conservation of 
mass and linear momentum. One might hope that all such descriptions would accordingly 
exhibit key energetic properties (nonlinear, exponential-type dissipation in the absence 
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of forcing and long-term stability under conditions of time dependent loading) consistant 
with the principle of energy conservation. These properties are intrinsic to real flows and 
the conventional, Eulerian Navier-Stokes equations. 

A completely general reference description of an incompressible, Newtonian fluid, which 
reconciles the differences between the so-called arbitrary Lagrangian Eulerian (A.L.E.) 
formulation of HUGHES, LiU and ZIMMERMAN |Q (deformation gradients absent) and 
that of SOULAIMANI, FORTIN, Dhatt and OUELLET @] (deformation gradients present, 
but use is problematic), is derived for the purposes of this investigation. The implications 
of the resulting description are investigated in the context of energy conservation in a 
similar, but broader, approach to that taken by others (eg. SiMO and Armero [H) for 
the conventional, Eulerian Navier-Stokes equations. 

The work subsequently focusses on establishing a class of time discretisations which 
inherit these self-same energetic properties irrespective of the time increment employed. 
The findings of this analysis have profound consequences for the use of certain classes 
of difference schemes in the context of deforming references. It is significant that many 
algorithms presently in use do not automatically inherit the fundamental qualitative 
features of the dynamics. An "updated" approach as a means of avoiding ever burgeoning 
deformation gradients which arise from the accumulated step-wise deformation of meshes 
and a still further simplified implementation are further topics explored. 

The main conclusions of this work are based on a new inequality and a number of lemmas. 
These lemmas are mainly concerned with the new convective term. The new inequality is 
used in place of where the Poincare-Friedrichs inequality might otherwise have limited the 
analysis. This analysis is extended in that non-zero boundaries, so-called free boundaries 
and time-dependent loads are considered. 



2 A Completely General Reference 

The implementation of most numerical time integration schemes would be problematic 
were a conventional Eulerian^ description of fluid motion to be used in instances involving 
deforming domains. The reason is that most numerical time integration schemes require 
successive function evaluation at fixed spatial locations. On the other hand meshes 
rapidly snarl when purely LagrangianQ descriptions are used. 

Eulerian and Lagrangian references are just two, specific examples of an unlimited number 
of configurations over which to define fields used to describe the dynamics of deforming 
continua. They are both special cases of a more general reference description, a de- 
scription in which the referential configuration is deformed at will. A deforming finite 
element mesh would be a good example of just such a deforming reference in practice. 
The transformation to the completely general reference involves coordinates where used 
as spatial variables only and the resultant description is therefore inertial in the same 

-'^ Eulerian or spatial descriptions are in terms of fields defined over the current configuration. 
^Lagrangian or material descriptions are made in terms of fields defined over a reference (material) 
configuration. 
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way as Lagrangian descriptions are. 
2.1 Notation 

Consider a material body which occupies a domain Q at time t. The material domain, 
Qq, is that corresponding to time t = to (the reference time, to, is conventionally, but 
not always, zero). A third configuration, Cl, which is chosen arbitrarily is also defined for 
the purposes of this work. The three domains are related in the sense that points in one 
domain may be obtained as one-to-one invertible maps from points in another. 




{Q. is spatially distorting Q., a 
deforming mesh, for example.) 

Figure 1: Schematic Diagram of Domains and Mappings Used in a Completely General 
Reference Description 

For any general function f{x,t), a function, f{x,t) = f{X*{x,t),t), can be defined 
in terms of the domains and one-to-one, invertible mappings illustrated in Figure 0. 
Similarly, fo{xo, t) = f{X{xo, t), t) can be defined. This notation can be generalised for 
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the component-wise definition of liiglier order tensors. Tlie key to understanding mucli 
of this work hes possibly in adopting a component-wise defined notation. 

In contrast to the function notation just estabhshed, the definition of the operators V 
and div is not based on V and div. They are instead the referential counterparts, that is 

~ d . d d d 

V = 7— and div = 7— — h 7— — h 



dx dxi dx2 dx 



3 



The notation A : B is used to denote the matrix inner product AijBij throughout this 
work, ( ■ , ■ )^2( . ) denotes the inner product and || ■ \\^2^ . ) the norm. 



2.2 Some General Results for Functions Defined on the Three 
Domains 

Three important results are necessary for the derivation of the completely general refer- 
ence description and these are presented below. 



The Material Derivative in Terms of a Completely General Reference 



The material derivative of any vector field v in terms of a completely general, reference 
is 



dv 



F ^v-v'^^] 



(1) 



where v^'^-^ is the velocity of the reference deformation, and F is the deformation gradient 
given by 

This result (taken from Hughes, Liu and Zimmerman g]) is obtained by recalling that 
the material derivative (total derivative) is the derivative with respect to time in the 
material configuration. Thus 



Dt 



d_ 

di 



{vi{X{xo,t),t)} 



dvi ^ dvi dXj 
dt dxj dt 



(2) 



A more practical expression is needed for (the velocity as perceived in the distorting 

at 

reference). This can be obtained by considering 

Xk{xo, t) = Xl{X{xo, t),t) (see Figure I on page |) 

so that 



dXk 



dt 



dXl 



Xo fixed 



dt 



+ 



dXl dX, 



X fixed 
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or 



dt dxk \ dt 



Xo fixed 



Substituting this expression into equation 
obtained. 



dt 



X fixed j 



the desired, suitably practicable result is 



An Element of Area in Terms of a Distorting Reference 

The second important result can be recalled from general continuum mechanics. Consider 
an element of area, size dA^ with an outward unit normal n. Then 

ndA = F'^NJdA (3) 

where dA and N denote the respective analogous size and outward unit normal of this 
element of area in the referential configuration (capital "n" so as to remain consistent 
with the notation, since Ni ^ rii in this case) and J = det F. This result is demonstrated 
in most popular textbooks on continuum mechanics (eg. Lai, Rubin and Krempl 0). 



The Kinematic Result J'q = Jq div v 

The material derivative of the Jacobian j7o is given by the relation 

Jo = Jo div V 

where j7o is defined as follows, 

Jo = det 



dx^]' 

This result is demonstrated in most popular textbooks on continuum mechanics (eg. Lai, 
Rubin and Krempl lIBi). 



2.3 Derivation of the Completely General Equation 

One way in which to derive a completely general reference description of an incompress- 
ible, Newtonian fiuid is to start with the balance laws in global (integral) form, and 
to make the necessary substitutions in these integrals. The desired numerical imple- 
mentation (similar to the conventional Navier-Stokes one which has been thoroughly 
investigated and found to be stable) is then obtained. 



Conservation of Mass 



Let r2(t) be an arbitrary sub- volume of material. The principle of conservation of mass 
states that 

d f 

— / pd^l = (rate of change of mass with time = 0) 
dt Jn(t) 
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d f 

-r / PoJodViQ = (reformulating in terms of the material 
at Jfio 

configuration, Qq.) 

f d 

/ TT {PoJ^o} d^o = (since limits are not time dependent in 
Juo at 

the material configuration.) (4) 
/ (poi7o + Poi7o) dQo = (by the chain rule) 



/ (p + pdivv)dri = (using the kinematic result Jq = JqASnv ) 
Jn(t) 



f I dv' doc ■ \ ~ ~ 

/ P + PTrr-TT^ Jdfl = (reformulating in terms of the distorting 
Jn(t) \ oxj oxi ) 

referential configuration, Q{t).) 
^ (p + pVi) : F ^ J = (integrand must be zero since the volume 

was arbitrary.) 

Thus, for a material of constant, non-zero density, 

Vv : ^ * = since J ^ (mappings are one-to-one and invertible) . 

Notice also that equation (H) implies 

^{PoJo} = (5) 
since the volume was arbitrary and the integrand must therefore be zero. 

Conservation of Linear Momentum (and Mass) 

The principle of conservation of linear momentum for an arbitrary volume of material 
Q{t) with boundary r{t) states that 

4- / pvd^l = [ pbdn + / crndA (6) 

dt Jn{t) Jn(t) Jr(t) 

where p is density, b is the body force per unit mass, cr is the stress, n the outward unit 
normal to the boundary and v is the velocity. The term on the lefthand side can be 
rewritten as follows: 



cL f a f 

— / pvdQ, = — poVoJodViQ (Reformulating in terms of the material 
dt Jn{t) dt Juo 

configuration, Qq-) 

d 

TT {pqVoJ^o} di^o (Since limits are not time dependent in 

Qo ut 

the material configuration.) 
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Q(t) 



n{t) 



pvJdVt 



pvdQ (The second term above is zero as a 

consequence of equation (^.) 
(Reformulating in terms of the dist- 
orting referential configuration, Cl.) 

JdCl (Using result 
on page 

where i) denotes the material derivative of v. The surface integral becomes 

crndA = I aF NJdA (Reformulating in terms of a distorting 
r{t) Jr{t) 

reference using result (Qj on page ^.) 



n{t) 



( dv 



div {ori^ J}dVt (By the divergence theorem). 



Finally, the term involving body force becomes 



n{t) 



phdVt 



n{t) 



pbJdQ 



(Reformulating in terms of a distorting 
reference.). 



Substituting these expressions into (|^), remembering that the volume used in the argu- 
ment was arbitrary and that the entire integrand must therefore be zero, the conservation 
principles of linear momentum and mass may be written in primitive form as 



p(^ + ViF'\v-v^-^f)\j 



pbJ + div P 



(7) 



and 



Vv : F * 







where P is the Piola-Kirchoff stress tensor of the first kind, P = crF J. In terms of 



the constitutive relation, cr = —pi + 2fiD, for a Newtonian fluid 
P = I -pi + p 



VvF + S/vF ^ 



F-'j 



smce 



The derivation of a variational formulation is along similar lines as that for the Navier- 
Stokes equations (the purely Eulerian description). For a fluid of constant density, the 
variational formulation 

. dv ~ ~ 



p I w ■ —JdQ + p w ■ Vv 
In at Jn 



F \v-v^^^] 



jdn 



p fw- bJdn + f pVw : F 'jdn - 2p f b{w) : t>{v)Jdn 
Jn Jn Jn 



wPNdT 



(9) 
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qVv -.F'dQ = (10) 

n 

is obtained, where q and w are respectively the arbitrary pressure and velocity of the 
variational formulation. 



2.4 Reconciling the Different Schools of Thought 

The equations (|^) and (^) are the completely general reference description of an incom- 
pressible, Newtonian fluid. They reduce to the so-called A.L.E. equations of Hughes, 
Liu and Zimmerman |^ for an instant in which spatial and referential configurations 
coincide. These simplified equations should, however, not be implemented where the im- 
plementation requires evaluation about more than one point within each time step (see 
Section ^ for a further, in-depth explanation). Under such circumstances the equations 
of Hughes et al. are an arbitrary Lagrangian Eulerian (A.L.E.) description in the 
very true sense (this is not surprising considering the equations have their origins in the 
arbitrarily, either Lagrangian or Eulerian programmes of HiRT, Amsden and CoOK 0). 
This fact is further borne out in observing that key energetic properties, consistant with 
the principle of energy conservation, are not automatically inherited by the equations of 
Hughes et. al. in the context of more general references. 

The momentum equations of SOULAIMANI, Fortin, Dhatt and Ouellet are flawed 

as a result of the mistaken belief that aF J is the Piola-Kirchoff stress tensor of the 
first kind (pg. 268 of Soulaimani et al.). Yet another problem is illustrated by 
rewriting the conventional incompressibility condition using the chain rule. The new 
incompressibility condition which arises is most certainly 

dvi dxj dvi dxi 

^TT^ = and not tt^t^ = 0. 

OXj OXi OXj OXj 

Further errors arising (eg. J omitted in the first term on the right hand side of the 
momentum equation, equation (10) on pg. 268 of Soulaimani et al.) make the use of 
these equations problematic. 



3 The Energetic Implications of a Deforming Refer- 
ence 



The effect of quantities parameterising reference deformation on key energetic properties 
- nonlinear, exponential-type dissipation in the absence of forcing and long-term stabil- 
ity under conditions of time dependent loading - is investigated in this section. These 
properties, 

K{v) < K{v ItJ e"^'^^* and lim sup K{v) < 



2Z/2C72 
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respectively (where K = ^p\\v\fj^2(^Q^ is the total kinetic energy), are intrinsic to real 



flows and the conventional, Eulerian Navier-Stokes equations (see Temam |]TD[, | |11| |, 
CONSTANTIN and FoiAS in and SiMO and Armero ^ in this regard). The effect of 
v'^'^-^ on the afore mentioned aspects of conservation of the quantity 



1 



Vj2 



L2(Q) 



is essentially what is being investigated, with a view to establishing a set of conditions 
under which the discrete approximation can reasonably be expected to inherit these self- 
same energetic properties. 

One might anticipate key energetic properties to be manifest only in instances involving 
a fixed contributing mass of material, whether its boundaries be dynamic, or not. An 
analysis of this nature only makes sense in the context of a constant volume of fluid 
which, for simplicity, will have material limits. 

Inequalities of the Poincare-Friedrichs type are a key feature of any stability analysis of 
this nature. Gradient containing terms need to be re-expressed in terms of energy. 
In the case of a "no slip" {v = 0) condition on the entire boundary the situation is 
straightforward, in that it is possible to use the standard Poincare-Friedrichs inequality: 
there exists a constant Ci > such that 

< Ci\\Vv\\l2 for all v e [H^i^T- 

The use of the classical Poincare-Friedrichs inequality is otherwise identified as a ma- 
jor limitation, even in the conventional Navier-Stokes related analyses. The Poincare- 
Friedrichs inequality is only applicable in very limited instances where the value for the 
entire boundary is stipulated to be identically zero. For boundary conditions of a more 
general nature, such as those encountered in this study, in which parts of the boundary 
may be either a free surface or subject to traction conditions, a more suitable inequality 
is required (notice that subtracting a boundary velocity and analising the resulting equa- 
tion is not feasible as the equations are nonlinear). The Poincare-Friedrichs inequality 
does, furthermore, not hold on subdomains of the domain in question and the constant 
is not optimal. 

Further investigation (communication 0) reveals a similar result, the so-called Poincare- 
Morrey inequality, holds providing the function attains a value of zero somewhere on the 
boundary. The Poincare-Morrey inequality states that a constant C2 > exists such that 

\\v\\l^ < C2\\Vv\\l2 for all v G [Hl^i^T- 

The proof of the Poincare-Morrey inequality is, however, similar to that of one of Korn's 
inequalities (see, for example, KiKUCHi and Oden [^]). In particluar, it is non-constructive, 
by contradiction and the constant cannot therefore be determined as part of the proof. 
Viewed in this light the forthcoming inequality amounts to a specification of the hypo- 
thetical constant in the Poincare-Morrey inequality for domains of a particular geometry. 
The particular types of geometry considered are those that arise in problems involving 
the motion of rigid bodies such as pebbles on the sea bed; thus a free surface is present, 
and the domain may be multiply connected. 
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Inequality 1 (A New "Poincare" Inequality) Suppose v is continuous and dif- 
ferentiable to first order and that v attains a maximum absolute value, c, on an included, 
finite neighbourhood of minimum radius Rmm about a point a;°'''g™ ( as depicted in Figure 




neighbourhood 



Figure 2: A Finite Neighbourhood of Minimum Radius /2min About a Point x 



origin 



If Q is a bounded, star-shaped (about a point a:;°"sinj|T| (domain in R'^, th 



en 



I'"IIl2(q) ^ 



{Rmax ^min) {Rjnax -^miny 



3 -Rmax Rm'm 

where i?max is the distance to the farthest point in Q from a:;°'''g™ 



Proof: Consider the change to spherical coordinates 

Vi{r, 6, (j)) = Vi{r sin 9 cos (j) — r sin Osmcj) — x™'^", r cos 6 — 2:3"^") 

centred on 0;°'''^™. Suppose the radial limits of the domain and neighbourhood are denoted 
Rh{0, (p) and Ra{0, 4>) respectively. By the fundamental theorem of integral calculus 



(yi{r,9,(j)) - u 



i \Ra{ 



lRa{e,4>) Or ^^'^'^"^^^ 



Ra{e,<i>) 4 or , 



l-r I 



' dvi 



^by which is meant that every point in the domain can be reached by a straight hne from x°^^^"^ that 
does not pass outside of 
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(by Schwarz inequality) 

/•i?max 1 rRb{d,<P) ( dVi V n 



max-' "rnin 



where V,ie,4>)= / C'dC 

jRa(e,4>) \dr J 

Integrating this result over that part of outside the neighbourhood (angular extent 
being e„(0) <e< 96(0) and $„ < < ^b) 

r^b fObW rRb{e,4>) , . 2 „ 

/ / / [vi{r,e,(j)) - Vi\R^(e^)] r sin 9drd9d(f) 



-'•'max -'•'min 

< ^ 1 p ^ / / V,{e, / r^dr sin ^cf^# 

^max-'t^min 



< 





3-Rmax-^^min 




(-^max 


~ -^min)(-Rmax 






3 Rmax Rmin 




(-^max 


~ Rimn)iRmax 


d3 \ 
-'''min/ 


3 Rmax Rmin 



' dvi \^ I ( dvA"^ 



1 / 



Or I \ do 

2" 

sin 9drd9d(j) 



sin 9 \ d4> ^ 

{Wvi) ■ {Wvi) sin 9drd9d(t). 



Changing back to the original rectangular coordinates and defining v |bndry to be a radially 
constant function throughout ^2 which takes the values of v |k„(6i,</,) for r = Ra{9, 0), 

/ {v,{x)-V,M'dn < (^■na.-^min)(g;La.-^^in) f (V^,(^)) . (V^,(^)) rfQ 

JQt <J-'''max-'''min Jfi* 



where is fl excluding the neighbourhood. Summing over i, 

I (t; - t; Ibndry) • (t; - Ibndry) rfl^ < (^"'^ " (^^a. " ^Lin) j 

Making use of either the Cauchy-Schwarz or triangle inequality, 

1^ , (-^max ^ -^min)(-Rmax ~ -^mi 



\V\ 



_ II I II \ ^ l.-n-max -'''minA-''"max -"-min/ ||Y7„,||2 

^J-^max-^min 
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and remembering that sup | v |i?„(e,(/,)|< c, 



L2(n*) 



< 



^- ) 

mm/ 



3-Rmax^min 



Consider the terms ||i^||j;^2 and ||c 



Il2- 



Comparing these terms under circumstances of 
sup I I < c leads to the conclusion that the inequality holds over the neighbourhood and 
that the inequality is therefore unaffected when the domain of integration is extended 
to include the neighbourhood. Of course, the radial extension of v |bndry can be used in 
place of c in instances where inclusion of the neighbourhood is not required. 

This inequality is similar to the Poincare-Friedrichs inequality when c = 0, but is ex- 
tended to a geometrical subclass of domains which have free and partly non-zero bound- 
aries. It has a further advantage in that the constant is an order of magnitude more 
optimal when used under the "no slip" Poincare-Friedrichs condition (under such con- 
ditions the domain can always be deconstructed into a number of subdomains in which 
-Rmin = |-Rmax)- The Poiucare-Friedrichs inequality is a special case of the above inequal- 
ity. The necessary lemma (below) follows naturally from the above inequality. 



Lemma 1 (Deviatoric Stress Term Energy) The kinetic energy satisfies the bound 



-K{v) 
P 



< 



^,^^^y where C is related to the constant in Inequality^, C > 0. 



Proof: If, in particular, v |bndry= in Inequality ||, 

llVvl 



\v\ 



L2(f7) 



< 



L2(n) 



C 



C- \ \v 



L2(n) 



< WDiv) 



lL2(r2) 



(The relationship between D and Vf arises in the context of the original equations 
involving diver. It is because 



= it;-- 

2 ^« 



(changing the order of differentiation) 
(divv = by incompressibility) , 



assuming, of course, that v is continuous and differentiable to first order.) Rewriting in 
terms of ^2 



-k(v) < 



D(v)j2 



L2(f7) 



The following lemma is vital to the deforming reference analysis in particular. It will 
form the basis to the next lemma and another (on page |l^) concerned with the time 
discrete analysis. 
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Lemma 2 (Basic to Lemmas |3| and |5]) The relation 

ii, (Vv)F'^wj\ = - (v, (\/u)F'^wj\ - (u (\/w : F^*) , vj) , 
is valid for 

weW = ^w : w = OorF iV ■ = on f| . 

Proof: Consider ii ■ {Vv)F ^wJ: 

UiVijF~i}wkJ = -UijViF'^^WkJ - UiVi{F~f}wkJ),j + {uiViF~^WkJ),j 

by the product rule. In the terms arising from {Fjf}wkJ)j, both Fjf}j and J,jF^}} van- 
ish under the condtions specified (in section for equations of Hughes, Liu and 
Zimmerman |^ to be a completely general reference description. Thus 

^iVijFji}wkJ = -UijViFff}wkJ - UiViFfi}wkjJ + {uiViFji}wkJ)j. 

Integrating over the domain Cl and applying the divergence theorem, 

u,{Vv)F'^wj) = -(v,{\/u)F'^wj) 



L2(Q) 

+ {u,v{f-'n-w)j)^^^^^ (11) 

The condition of this lemma dictates the manner in which the reference must deform 
to ensure that the equation will inherit the desired energetic properties. This lemma is 
crucial to the deforming reference analysis. The lemma immediately below will facilitate 
the elimination of the convective energy rate in the forthcoming analysis. 

Lemma 3 (Convective Energy Rate) The relation 

p(v,{Vv)F'Uv-VAj) _ = -l-p(v,v—\ 



I L2(n) 2 \ ' dt / 



L2(n) 



is valid in instances where a purely Lagrangian description is used to track free boundaries 
and/or boundaries are of a fixed impermeable type. 

Proof: In instances where a purely Lagrangian description is used to track free bound- 
aries, v — v^'^-' vanishes, as does F N -v at fixed impermeable boundaries. The condition 
at the boundary for Lemma ^ is therefore satisfied. Thus the term 

- p(v,iVv)F'^ (v-v'-'f) j) _ = p/v,{Vi)F'\v-v'^'^) j) 



+p (t) (V {v - v'^'f) : ^ *) > "i^J 



L2(n) 
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\ ^ ^ / L2{Q) 

= -ip(*(Vi"/:F-'),*j)^.^^^ 

1 

C?t/ ~f — ~f — * 

since --- = Jdivv'^'^^ (which is J'Vv'^'^-' : F ) in the same vein as j7o = j7o div v (the 
kinematic result on used earher). 

This lemma concludes the preliminaries required for the deforming reference energy anal- 
ysis. 



3.1 Exponential Dissipation in the Absence of Forcing 

The issue of whether nonlinear, exponential-type dissipation in the absence of forcing is 
a property intrinsic to the deforming reference description is resolved as follows. 



Theorem 1 (Exponential Dissipation in the Absence of Forcing) A sufficient 
condition for the completely general reference description to inherit nonlinear, exponential 
type energy dissipation 

k{v) < k{v lij e-^^^^* 

1 . ~i 2 

(where K{v) = -p vJ^ ^2 ^ ^ absence of forcing (an intrinsic feature of real flows 

and the conventional, Eulerian Navier-Stokes equations) is that the reference moves in a 
purely Lagrangian fashion at free boundaries. 



Proof: The first step towards formulating an expression involving the kinetic energy is 
to substitute v for w in the variational momentum equation (H) on page |^. Then 

p(v,^j\ = (pVv,F''j)^.-2p(D(v).b(v)j)^. 

-p /v, {Vv)F'^ (i; - v'-'f) j\ ^ _ 



+p(v,bJ) - + (v,PN) , - . (12) 



The term containing the pressure, that is 
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vanishes as a result of incompressibility (equation (H)). The order of integration and 
differentiation are interchangeable (limits are time-independent in the reference which 
tracks the free boundary perfectly - a description which becomes fully Lagrangian at 
boundaries was stipulated). Equation (|l2|) can be rewritten 



d 
di 



Vj'2 



2 I dj\ 







) = -2// 


d3 



2 



+p(v,bJ) , - + (i,PN] 

as a result. The conditions of Lemma |^ are also satisfied for a description which becomes 
fully Lagrangian at free boundaries and an expression 

dk{v) 



dt 



= -2/i DJ2 
is therefore obtained, where K = -p 



, - +p{v,hJ) ^ - + (v,PN) ^ - 



- is the total kinetic energy. Using Lemma 



< -2uCk(v)+p(v,bj) ^ . + (v,Pn] _ 



(13) 



Equation (|T3D has a solution of the form 

K < k{v 



-2vCt 



in the absence of forcing ( "no forcing" =^ b = PN = 0), providing a purely Lagrangian 
description is used at free boundaries. 

A nonlinear, exponential-type energy dissipation in the absence of forcing is therefore an 
intrinsic property of the completely general reference description. This contractive flow 
property is also an intrinsic property of the conventional Navier-Stokes equations. 



3.2 Long— Term Stability under Conditions of Time— Dependent 
Loading 

The formulation of suitable load and free surface bounds is necessary before the issue 
of long-term stability (L^-stability) under conditions of time-dependent loading can be 
resolved. The following lemma facilitates the formulation of load and free surface bounds. 



Lemma 4 (Force, Free Surface Bounds) The inequality 



p{v,bJ) , +(v,PN) , < 



p 



Vj2 



I 1 1 ~ I |2 



2uC 



p 



bj2 



L2(n) 



+ 



L2(r) 



PN 



L2(f) 



holds where uC is a constant, uC > 0. 
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Proof: In terms of the Cauchy-Schwarz inequality, 



v,bJ) , - < 

L2(n) - 

< 



~1 




6J5 









L2(n) 



+ 



bj2 



by Young's inequality. Similarly, 



L2(r) 



< \\v\\r2m^ + 



PAT 



L2(n) 



L2(f) 



for //C > 



for uC > 0. 



This done, the mathematical machinery necessary to the long-term stability analysis is 
in place. 



Theorem 2 (Long-Term Stability) A sufficient condition for the completely gen- 
eral reference description to inherit the property of long-term stability 

limsupir(i;)<^^ 

under conditions of time-dependent loading (an intrinsic feature of real flows and the 
Navier-Stokes equations), where this time-dependent loading and the speed of the free 
surface is bounded in such a way that 



bJ^ 



PN 



L2(f) 



+ iy'C'\\v 



iL2(r) 



is that the description becomes purely Lagrangian at free boundaries. 
Proof: Using Lemma ^ in equation (|13|), then applying the above bound. 



dk{v) ^-^^^ ^ 

Using the Gronwall lemma (see HiRSCH and Smale 0) leads to the differential inequality 

M2 



dK{v) .y. ct 



dt 



2uC 



which, when solved, yields 



K{v) < e-'''''K{v\t=to) + il-e 



-vCt 



M2 



2l^2C2 ' 



This in turn implies 



M2 

limsupK(v) < 7-^7^. 



The preceding analyses lead to natural notions of nonlinear dissipation in the absence 
of forcing and long-term stability under conditions of time-dependent loading for the 
analytic problem. 
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The Energetic Implications of the Time Discreti- 
sation 



This section is concerned with estabhshing a class of time discretisations which inherit 
the self-same energetic properties (nonlinear dissipation in the absence of forcing and 
long-term stability under conditions of time dependent loading) as the analytic problem, 
irrespective of the time increment employed. In this section a generalised, Euler difference 
time-stepping scheme for the completely general reference equation is formulated and 
the energetic implications investigated in a similar vein as the analytic equations in the 
previous section. 

This stability analysis is inspired by the approach of others to schemes for the conventional 
Navier-Stokes equations. The desirability of the attributes identified as key energetic 
properties is recognised and they have been used as a benchmark in the analysis of 
various of the conventional, Eulerian Navier-Stokes schemes by a host of authors. Related 
work on the conventional, Eulerian Navier-Stokes equations can be found in a variety of 



references, for example Temam [11| and SiMO and Armero 



The analyses presented here are extended, not only in the sense that they deal with the 
completely general reference equation, but also in that non-zero boundaries, so-called 
free boundaries and time-dependent loads are able to be taken into account (the former 



two as a consequence of the new inequality on page |I0|) . The findings of this work have 
profound consequences for the implementation of the deforming reference equations. It 
is significant that many algorithms used for long-term simulation do not automatically 
inherit the fundamental qualitative features of the dynamics. 



A Generalised Time— Stepping Scheme 

An expression for a generalised Euler difference time-stepping scheme can be formulated 
by introducing an "intermediate" velocity 

Vn+a = a V |t+At +(1 " ") ^> |t for a e [0, 1] (14) 

to the variational momentum equation (equation (y) on page |^ where v \t and v \t+At 
are the solutions at times t and t + At respectively. At being the time step. It is in this 
way that a generalised time-discrete approximation of the momentum equation 

+p (w, b 

\ I L-'{iln + a) \ I (Vn + a) 

is derived, where { ■ ) L2(p^ ^ denotes the inner product over the deforming domain 

at time t + a At. r„+a, Fn+a, Jn+a, Dn+a, Pn+a, and bn+a are likewise defined to be 
the relevant quantities evaluated at time t + aAt. 
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It will presently become clear that it makes sense to perform the analyses for the time- 
discrete equation in the context of divergence free rates of reference deformation only. 
This is since relevant energy terms are not readilly recovered from the time-discrete 
equations for deforming references in general. This investigation is accordingly restricted 
to a subclass of reference deformations in which "reference volume" is conserved. This 
is for reasons of expedience alone and the subclass of deformations is thought to be 
representative. 



Assumption 1 The assumptions Jn = Jn+a and Jn+i = Jn+a are made so that 



K{vn) = -p 



~i 


2 


1 


_ 1 


VnJn 













and 



K{vn+i) = -p 



L2{Q„+i) 2 

(by equation ([T^j and since the volume of material over which integration is being per- 
formed is constant). 



J J d J 

Remark: Notice that = Jdivi;!^'^^, the discrete form of — — = Jdivi;'"'^-'^, can 

At ^ dt 

consequently be rewritten as 

div vl,%„ = 



under the conditions of the above assumption. It is for the practical expedience afforded 
by Assumption |l| alone that this analysis is limited to instances in which divu^^^^ = 0. 

The following lemma will establish that the rate of energy change associated with the 
convective term vanishes as a result of the assumption. 



Lemma 5 (Discrete Convective Energy Rate) The discrete convective term 

-p(w, {VVn+a)F^l^ (Vn+a " K+a) Jn+a) , 

vanishes under circumstances of div f^+^j = and a purely Lagrangian description is 
used at free boundaries (alternatively boundaries are of the fixed, impermeable type). 



Proof: The operator ( ■ , (V 



)F ^wj\ 



w eW = {w : (Vw) : F 



on Cl; w 



is skew-symmetric for 



or F ^A/" ■ -0; = on r} 



by equation ( pA]) on page 13. In instances where a purely Lagrangian description is 
used to track free boundaries ■&„+q, — v^la vanishes. At fixed, impermeable boundaries 
vanishes. The condition at the boundary is therefore satisfied, under 



F * N ^ ■ V 
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all of the afore-mentioned circumstances. Apply the stipulated condition Vv'^'^-' : F =0 
and set w = Vn+a — ^n+a etc. 

Remark: Recall that in the investigation of the analytic problem, a term arising from the 
manipulation of the acceleration containing term (the term containing the rate of change 
of the Jacobian) cancelled with the convective energy. It is therefore not surprising that 
assumptions pertaining to the acceleration containing term (in particular to the rate of 
change of the Jacobian) in the discrete problem will, once made, also be necessary for 
the corresponding discrete convective energy term to vanish (reffering to the div v^^^ = 
condition of Lemma ^). This is a good prognosis for the energetic behaviour of the 
discrete problem in circumstances of reference deformations excluded by Assumption |l]. 

This concludes the preliminaries required for the analysis of the time-discrete equation. 



4.1 Nonlinear Dissipation in the Absence of Forcing 

The following analysis establishes a class of time-stepping schemes which exhibit nonlin- 
ear dissipation in the absence of forcing regardless of the time increment employed. 

Theorem 3 (Nonlinear Dissipation in the Absence of Forcing) Suppose that 
the description is pure Lagrangian at any free boundaries and that the deformation rate 
of the reference is divergence free. A sufficient condition for the kinetic energy associated 
with the generalised class of time-stepping schemes to decay nonlinearly 



k{vn+i) - Kivn) < -At2/i 



in the absence of forcing and irrespective of the time increment employed, is that the 
scheme is as, or more, implicit than central difference. That is 

1 



Proof: Expressing the intermediate velocities v^^i and i7„+o in terms of equation ( [I^ ) 
and subtracting, the result 

is obtained. The first step towards formulating an expression involving the kinetic energy 
of the generalised time stepping-scheme (^) is to replace the arbitrary vector, w, with 
Vn+a- By further substituting (|16D into ( plSf ) and eliminating the pressure containing term 



in a similar manner to that in Theorem |l|, an expression involving the difference in kinetic 
energy over the duration of a single time step is obtained. 



The vector v — v"^^^ vanishes in instances where a purely Lagrangian description is used 
to track free boundaries. The quantity F^^^Nn+a ■ Vn+a vanishes where boundary con- 
ditions are of a fixed impermeable type. The condition at the boundary for Lemma ^ is 
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therefore satisfied. Incompressibility and a restriction on reference deformations to those 
for which divvl^J^^ is zero ensure that the remaining Lemma ^ condition is satisfied. 



The equation 



r2 



-At 2/i 



2 



is then obtained. Since it is assumed that there is no forcing 

1 



(17) 



K{vn+i) - K{vn) < -pi a 



-At 2/i 



Vn+l - Vn) J,' 



D(v 



n+a) Jn+a 



Thus the kinetic energy inherent to the algorithmic fiow decreases nonhnearly in the 
absence of forcing, irrespective of the time increment employed and for arbitrary initial 
conditions provided that 

1 



a > 



and 



div vl^l^ 



0. 



The former requirement translates directly into one specifying the use of schemes as, 
or more, implicit than central difference. Only for descriptions which become fully La- 
grangian at free boundaries can it be guaranteed that energy will not be artificially 
introduced by way of the reference. 

Remark: Notice (by Lemma P that for a = | an identical rate of energy decay 

K{Vn+l) - K{Vr 



At 



< 



is obtained for the discrete approximation as was obtained for the equations. 



4.2 Long— Term Stability under Conditions of Time— Dependent 
Loading 

This second part of the time-discrete analysis establishes a class of time stepping schemes 
which exhibit long-term stability under conditions of time dependent loading irrespective 
of the time increment employed. The following lemma is necessary to the analysis and 
is concerned with devising a bound for the energy at an intermediate point in terms of 
energy values at either end of the time step. 

Lemma 6 (Intermediate Point Energy) The following bound applies 

k{vn+a) >a{a-c + ac) k{vn+i) + (1 - a) (^1 - a - k{vn) 
where c is some constant, c > 0. 
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Proof: By Young's inequality 



11 1^ 



< 



+ 



2 



for c > 0. Writing K{vn+a) explicitly in terms of the "intermediate" velocity definition, 
leads to 



> a^kivn+i) + (1 - afkivn) 
-2q;(1 - a) v^+iJn+a 

> a [a — (1 — a)c] K{Vn+i) + (1 — a) (1 — a) — — K{Vr 



V P 



a 



using equation ([TsD. The optimal choice of the constant c is established farther on. 

The following theorem establishes a class of time-stepping schemes which exhibit long- 
term stability under conditions of time-dependent loading regardless of the time incre- 
ment employed. 



Theorem 4 (Long-Term Stability) Suppose that the description is pure Lagrangian 
at any free boundaries and that the rate at which the reference is deformed is divergence 
free. A sufficient condition for the algorithmic flow to exhibit long-term stability under 
conditions of time-dependent loading (intrinsic to real flows and the Navier-Stokes equa- 
tions), assuming this time-dependent loading and the speed of the free surface is bounded 
in such a way that 





_1 


2 




p 




+ 


P ]\f 

n+a-^ ^ n+a 











n+a\\L2{r„+^) 



IS 



1 

a > -. 
2 



Proof: Substituting Lemma ^ (page |T5D and Lemma |I| (page [T^) into equation ([T 
applying the above bound and choosing a > | one obtains 



At 



< 



2^' 



From this point on the argument used is identical to that of SiMO and Armero |^ for 
the conventional, Eulerian Navier-Stokes equations. Substitution of Lemma || leads to a 
recurrence relation, 



T<^(~ ^ <, l-^C(l-a)(l-a-f)At g,^. 

1 + uCa[a — c + ac)At 



M^At 



2uC [1 + uCa{a - 1 + ac) At] ' 
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Using this recurrence relation to take cognisance of the energy over all time steps, 
'1 - uC{l - a){l - a - ^)At' 



+1J 



< 



1 + vCai^a — c + ac) At 



n-l 

E 



2vC [1 + pCa{a - c + ac) At] 



;i-i/C(l-a)(l-a-f)At) 
1 + uCa{a — c + ac) At 



(19) 



is obtained. An infinite geometric series which converges so that 

M^At 



lim sup iir(v„+i) < 



2uC [1 + uCa{a - c + ac) At] 



1 



(l-i/C(l-a)(l-a-f)At) 
1 + vCa^a — c + ac) At 

M2 



n -1 



2yC [uCa{a -c + ac) + iyC{l - a)(l - a - f ) 
results, providing the absolute ratio of the series is less than unity. That is 



1 + i/Ca(Q; — c + ac) At 



< 1. 



Therefore either 



a 



or 



-1 - uCa{a - c + ac) At < 1 - uC{l - a) ( 1 - a ) At 



1 - iyC{l - a) [l-a-^] At <1 + uCa{a -c + ac)At 



(20) 



in order for the bound to exist. Notice, furthermore, that for this desired convergence to 
be unconditional (regardless of the time increment employed) requires 

a-c + ac>0. (21) 

The denominator in the series ratio might otherwise vanish for some value of At. 



For a G 



2' 



equation ([20| ) and equation ( pTj) together imply 
(I -a) 



<c< 



a 



a 



:i-a) 



which in its turn implies 



:i-a) 



a 



< 



a 



[l-aY 
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The choice of the parameter a > ^ therefore leads to an infinite geometric series which 
forms the desired upper bound. The minimum value of this bound occurs for c chosen 
according to 

inf 



(i^<,<^ uCa(a -c + ac) + uCil - a){l - a - ^) pC{2a - If ' 

The value of this upper bound, which occurs for the choice of the parameter a > ^, is 
then 

lim sup ii'(t)„+i) < 



2z^2C2(2a- 1)2- 

In this way one arrives at a class of algorithms which are unconditionally (irrespective of 
the time increment employed) stable. 

Remark: Notice that for a = 1 one obtains an identical energy bound for the discrete 
approximation as was obtained for the equations. 



5 An "Updated" Approach and a Simplified Imple- 
mentation 

Ever burgeoning deformation gradients accumulate for a straight forward implementa- 
tion of the equations. Using an "updated" approach is one way of coping with this 
otherwise rather daunting prospect. An "updated" approach is the result of a little, 
well-worthwhile lateral thinking. An "updated" approach amounts to choosing a new 
referential configuration during each time step. 

In the case of time stepping schemes based about a single instant (eg. the generalised 
class of Euler difference schemes investigated in Section ^) a considerably simplified im- 
plementation can further be achieved by a particularly appropriate choice of configura- 
tions. Making the choice of a referential configuration which coincides with the spatial 
configuration at the instant about which the time stepping scheme is based allows the 
deformation gradient to be omitted altogether (the deformation gradient is identity under 
such circumstances). For such implementations (those which require evaluation about a 
single point only) no error arises from the use of the equations cited in Hughes, Liu 
and Zimmerman |||, 

p(^^ + Vv{v-v'''f)^ = p6 + diver (22) 
dwv = 0. (23) 

These equations are not valid for any, arbitrary choice of reference or if the implemen- 
tation requires the equation to be evaluated at more than one point within each time 
step (eg. a Runge-Kutta or finite-element-in-time scheme). It is important to remember 
that in a discrete context the reference configuration is fixed for the duration of the entire 
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time increment. Although the referential configuration is hypothetical and can be chosen 
arbitrarily for each time step, once chosen it is static for the duration of the entire time 
step. Once the coincidence of configurations is ordained at a given instant, F is defined 
by the deformation, both before and after, and must be consistant. 

There would seem to be no reason why one would wish to define the deformation about 
a configuration other than that at the instant about which the implementation is based 
(assuming the implementation used is indeed based about a single point eg. a finite 
difference) thereby involving deformation gradients. Resolving the resulting difficulties 
associated with the deformation gradients by means of a perturbation seems unnecessarily 
complicated in the light of the above reasoning. 



6 Conclusions 



The correct equations, which describe the motion of an incompressible, Newtonian fiuid 
and which are valid for a completely general range of reference deformations, are equations 
(0) and (P). For implementations requiring the equations to be evaluated about a single 
instant within each time step only (eg. finite differences), the deformation gradients may 
be assumed identity i.e. the equations of Hughes, Liu and Zimmerman M (equations 



(H) and (ID) will suffice. 



In this work it is shown (as was hoped) that nonlinear, exponential-type dissipation in the 
absence of forcing and long-term stability under conditions of time dependent loading 
are properties automatically inherited by deforming reference descriptions. The single 
provisor is that such descriptions become fully Lagrangian at any moving boundaries. 
These properties are intrinsic to real flows and the conventional, Eulerian Navier-Stokes 
equations. 

Relevant energy terms are not readily recovered from the time-discrete equations for 
deforming references in general. Only for divergence free rates of reference deformation 
which become fully Lagrangian at free boundaries could it consequently be guaranteed 
that energy would not be artificially introduced to the algorithmic fiow by way of the 
reference. The divergence free assumption was made for reasons of expedience alone and 
the limitations of the time-discrete analysis are consequently not expected to detract from 
the use of the method in any way. This is especially so when it is considered that, a term 
arising from the manipulation of the acceleration containing term (the term containing the 
rate of change of the Jacobian) cancelled with the convective energy in the investigation 
of the analytic problem and that assumptions pertaining to the acceleration containing 
term (in particular to the rate of change of the Jacobian) in the discrete problem were, 
once made, also necessary for the corresponding discrete convective energy term to vanish 
(reffering to the div v*"^-^ = condition of Lemma |^). If one were to be overly cautious on 
this basis one would be faced with the additional challenge of enforcing a fully Lagrangian 
description for nodes situated on any free boundaries, while deforming elements would 
be required to deform at a rate which is divergence free. Such a totally divergence 
free description may, however, not be possible. An alternative strategy would be to use 
a fully Lagrangian description. Both the purely Lagrangian and purely Eulerian fluid 
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descriptions have divergence free rates of distortion. 

There are inherent problems with using certain classes of time-stepping schemes and 
the use of finite difference schemes more implicit than central difference is consequently 
advocated. Such differences exhibit the key energetic properties (nonlinear, exponential- 
type dissipation in the absence of forcing and long-term stability under conditions of time 
dependent loading) irrespective of the time increment employed. A backward difference 
is the obvious choice. Calculations at time t + aAt would require an intermediate mesh 
and associated quantities for instances in which a ^ 1 (since a > |). 

The author recommends a strategy in which a predominantly Eulerian description is 
used, where possible, for the bulk of the problem (from an efficiency point of view) and 
the completely general reference description for the remainder is appropriate. Purely 
Eulerian descriptions have the advantage of a "one off" finite element construction and 
involve none of the hazards of a badly distorted reference. 
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